!// Dmitry.Konovalov@jcu.edu.au Dmitry.A.Konovalov@gmail.com, 2013

module gamma_mod;
use dbg_mod;
implicit none;
private;
public DGAMMA, DLGAMA, test_gamma_mod;
contains

subroutine test_gamma_mod();
  character(LEN=64) :: name = "test_gamma_mod";
  call dbg_in(name); 
  call test_gamma();
  call dbg_out(name); 
end subroutine

!! from http://www.sourceware.org/bugzilla/show_bug.cgi?id=5159
!#include <math.h>
!#include <stdio.h>
!#define N_MAX 34
!#undef M_PI
!#define M_PI 3.14159265358979323846
!int main()
!{
!  float c[N_MAX+1];
!  float xs, ts;
!  float diff;
!  int i;
!
!  c[0] = 1.;
!  for (i=1; i<=N_MAX; i++)
!    c[i] = (2*i-1)*c[i-1]*0.5;
!
!  for (i=1; i<=N_MAX; i++)
!    {
!      xs = i + 0.5f;  /* argument to the gamma function */
!      ts = c[i] * sqrtf(M_PI);
!      diff = fabsf(tgammaf(xs)-ts)/ts;
!      if (diff > 1e-6)
!        printf("Arg: %g, Exact: %12.8g, tgamma: %12.8g, diff: %12.8g\n",
!              xs, ts, tgammaf(xs), diff);
!    }
!  return 0;
!}
SUBROUTINE test_gamma();  
  character(len=1024) :: mssg;
  integer, parameter :: N_MAX = 34;
!  float c[N_MAX+1];
  real(DP) :: c(0:N_MAX+1)
  real(DP) xs, ts, diff, di, logts;
  integer i;

  !call dbg_assert_abs_diff_LE("M_PI", PI, M_PI, EPS16);

! c[0] = 1.;
  c(0) = D1;
! for (i=1; i<=N_MAX; i++)
  do i = 1, N_MAX
    di = real(i, DP);
!   c[i] = (2 * i - 1) *c[i-1]*0.5;
    c(i) = (D2*di - D1)*c(i-1)*D05;
  end do

!  for (i=1; i<=N_MAX; i++)
  do i = 1, N_MAX
    di = real(i, DP);    !call dbg_dp("di", di);
!   xs = i + 0.5f;  /* argument to the gamma function */
    xs = di + D05; 
!   ts = c[i] * sqrtf(M_PI);
    ts = c(i) * SQRT(PI);
    
    !call dbg_dp("SQRT(M_PI) ", SQRT(M_PI));
    !call dbg_dp("DSQRT(M_PI)", DSQRT(M_PI));
          
    diff = ABS(dgamma(xs)-ts)/ts;
    write(mssg, *)  "Arg=", xs, ", Exact=", ts, ", dgamma=", dgamma(xs) &
    &, ", diff=", diff;
    call dbg_str(TRIM(mssg));
    call dbg_assert_abs_diff_LE(TRIM(mssg), diff, D0, EPS15);
    
    logts = LOG(ts);
    diff = ABS(DLGAMA(xs) - logts) / logts;
    write(mssg, *)  "Arg=", xs, ", Exact LOG(GAMMA)=", logts, ", dgamma=", DLGAMA(xs) &
    &, ", diff=", diff;
    call dbg_str(TRIM(mssg));
    call dbg_assert_abs_diff_LE(TRIM(mssg), diff, D0, EPS13);
    
   end do
end SUBROUTINE




!! original from gen.f
!      FUNCTION DGAMMA(X)
!C
!      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
!C
!      CHARACTER*(*) NAME
!      PARAMETER(NAME='GAMMA/DGAMMA')
!C
!      CHARACTER*80 ERRTXT
!
!      DIMENSION C(0:15)
!
!      DATA C( 0) /3.65738 77250 83382 44D0/
!      DATA C( 1) /1.95754 34566 61268 27D0/
!      DATA C( 2) /0.33829 71138 26160 39D0/
!      DATA C( 3) /0.04208 95127 65575 49D0/
!      DATA C( 4) /0.00428 76504 82129 09D0/
!      DATA C( 5) /0.00036 52121 69294 62D0/
!      DATA C( 6) /0.00002 74006 42226 42D0/
!      DATA C( 7) /0.00000 18124 02333 65D0/
!      DATA C( 8) /0.00000 01096 57758 66D0/
!      DATA C( 9) /0.00000 00059 87184 05D0/
!      DATA C(10) /0.00000 00003 07690 81D0/
!      DATA C(11) /0.00000 00000 14317 93D0/
!      DATA C(12) /0.00000 00000 00651 09D0/
!      DATA C(13) /0.00000 00000 00025 96D0/
!      DATA C(14) /0.00000 00000 00001 11D0/
!      DATA C(15) /0.00000 00000 00000 04D0/
!
!      U=X
!      IF(U .LE. 0) THEN
!       WRITE(ERRTXT,101) U
!       CALL MTLPRT(NAME,'C302.1',ERRTXT)
!       H=0
!       GO TO 9
!      ENDIF
!    8 F=1
!      IF(U .LT. 3) THEN
!       DO 1 I = 1,INT(4-U)
!       F=F/U
!    1  U=U+1
!      ELSE
!       DO 2 I = 1,INT(U-3)
!       U=U-1
!    2  F=F*U
!      END IF
!      H=U+U-7
!      ALFA=H+H
!      B1=0
!      B2=0
!      DO 3 I = 15,0,-1
!      B0=C(I)+ALFA*B1-B2
!      B2=B1
!    3 B1=B0
!    9 DGAMMA=F*(B0-H*B2)
!      RETURN
!  101 FORMAT('ARGUMENT IS NEGATIVE = ',1P,E15.1)
!      END 

      FUNCTION DGAMMA(X)
      !IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      integer I;
      real(DP) X, U, H, F, DGAMMA, ALFA, B1, B2, B0; 
      CHARACTER*(*) NAME
      PARAMETER(NAME='GAMMA/DGAMMA')
      CHARACTER*80 ERRTXT

      real(DP) :: C(0:15)
      DATA C( 0) /3.65738772508338244_DP/
      DATA C( 1) /1.95754345666126827_DP/
      DATA C( 2) /0.33829711382616039_DP/
      DATA C( 3) /0.04208951276557549_DP/
      DATA C( 4) /0.00428765048212909_DP/
      DATA C( 5) /0.00036521216929462_DP/
      DATA C( 6) /0.00002740064222642_DP/
      DATA C( 7) /0.00000181240233365_DP/
      DATA C( 8) /0.00000010965775866_DP/
      DATA C( 9) /0.00000000598718405_DP/
      DATA C(10) /0.00000000030769081_DP/
      DATA C(11) /0.00000000001431793_DP/
      DATA C(12) /0.00000000000065109_DP/
      DATA C(13) /0.00000000000002596_DP/
      DATA C(14) /0.00000000000000111_DP/
      DATA C(15) /0.00000000000000004_DP/

      U=X
      IF(U .LE. 0) THEN
        WRITE(ERRTXT,101) U
!       CALL MTLPRT(NAME,'C302.1',ERRTXT)
       call dbg_assert_true(ERRTXT, U > 0);
       H = D0;
       GO TO 9
      ENDIF
    8 F = D1;
      IF(U .LT. 3) THEN
        !call dbg_dp("   (4-U)=", (4-U));        
        !call dbg_int("INT(4-U)=", INT(4-U));
        DO I = 1, INT(4-U)
          F = F / U;
          U = U + D1;
        end do
      ELSE
        !call dbg_dp("   (U-3)=", (U-3));
        !call dbg_int("INT(U-3)=", INT(U-3));
        DO I = 1, INT(U-3)
         U = U - D1;
         F = F * U;
        end do
      END IF
      H = U + U - D7;
      ALFA = H + H;
      B1 = D0;
      B2 = D0;
      
      !call vec_dbg("C(1:15)=", C);
      DO I = 15,0,-1      
        B0 = C(I) + ALFA * B1 - B2; 
!        call dbg_dp("CI=", C(I)); 
!        call dbg_dp("B0=", B0);
        B2 = B1;
        B1 = B0;  
      end do
    9 DGAMMA = F * (B0 - H * B2);
      RETURN
  101 FORMAT('ARGUMENT IS NEGATIVE = ',1P,E15.1)
      END function




!      FUNCTION DLGAMA(X)
!C
!      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
!C
!      CHARACTER*(*) NAME
!      PARAMETER(NAME='ALGAMA/DLGAMA')
!C
!      DIMENSION P1(7),Q1(7),P2(7),Q2(7),P3(7),Q3(7),C(5)
!
!      PARAMETER (Z1 = 1, HF = Z1/2, HF1 = 1+HF)
!      CHARACTER*80 ERRTXT
!      DATA P1
!     1/+3.84287 36567 45991D+0, +5.27068 93753 00983D+1,
!     2 +5.55840 45723 51531D+1, -2.15135 13573 72570D+2,
!     3 -2.45872 61722 29242D+2, -5.75008 93603 04123D+1,
!     4 -2.33590 98949 51284D+0/
!      DATA Q1
!     1/+1.00000 00000 00000D+0, +3.37330 47907 07074D+1,
!     2 +1.93877 84034 37713D+2, +3.08829 54973 42428D+2,
!     3 +1.50068 39064 89095D+2, +2.01068 51344 33395D+1,
!     4 +4.57174 20282 50299D-1/
!      DATA P2
!     1/+4.87402 01396 83863 6D+0, +2.48845 25168 57407 6D+2,
!     2 +2.17973 66058 89591 5D+3, +3.79751 24011 52511 8D+3,
!     3 -1.97780 70769 84164 6D+3, -3.69298 34005 59128 2D+3,
!     4 -5.60177 73537 80387 7D+2/
!      DATA Q2
!     1/+1.00000 00000 00000 0D+0, +9.50999 17418 20893 8D+1,
!     2 +1.56120 45277 92863 5D+3, +7.23400 87928 94807 1D+3,
!     3 +1.04595 76594 05895 9D+4, +4.16994 15153 20023 1D+3,
!     4 +2.76785 83623 80410 1D+2/
!      DATA P3
!     1/-6.88062 40094 59425D+3, -4.30699 69819 57098D+5,
!     2 -4.75045 94653 43956D+6, -2.94234 45930 32234D+6,
!     3 +3.63218 04931 54257D+7, -3.35677 82814 54576D+6,
!     4 -2.48043 69488 28593D+7/
!      DATA Q3
!     1/+1.00000 00000 00000D+0, -1.42168 29839 65146D+3,
!     2 -1.55528 90280 85353D+5, -3.41525 17108 01107D+6,
!     3 -2.09696 23255 80444D+7, -3.45441 75093 34395D+7,
!     4 -9.16055 82863 71317D+6/
!      DATA C
!     1/ 1.12249 21356 561D-1,  7.95916 92961 204D-2,
!     1 -1.70877 94611 020D-3,  9.18938 53320 467D-1,
!     2  1.34699 05627 879D+0/
!
!
!      ENTRY DLOGAM(X)
!
!      IF(X .LE. 0) THEN
!       H=0
!       WRITE(ERRTXT,101) X
!       CALL MTLPRT(NAME,'C304.1',ERRTXT)
!      ELSE IF(X .EQ. 1 .OR. X .EQ. 2) THEN
!       H=0
!      ELSE IF(X .LE. HF) THEN
!       Y=X+1
!       AP=P1(1)
!       AQ=Q1(1)
!       DO 2 I = 2,7
!       AP=P1(I)+Y*AP
!    2  AQ=Q1(I)+Y*AQ
!       H=-LOG(X)+X*AP/AQ
!      ELSE IF(X .LE. HF1) THEN
!       AP=P1(1)
!       AQ=Q1(1)
!       DO 3 I = 2,7
!       AP=P1(I)+X*AP
!    3  AQ=Q1(I)+X*AQ
!       H=(X-1)*AP/AQ
!      ELSE IF(X .LE. 4) THEN
!       AP=P2(1)
!       AQ=Q2(1)
!       DO 4 I = 2,7
!       AP=P2(I)+X*AP
!    4  AQ=Q2(I)+X*AQ
!       H=(X-2)*AP/AQ
!      ELSE IF(X .LE. 12) THEN
!       AP=P3(1)
!       AQ=Q3(1)
!       DO 5 I = 2,7
!       AP=P3(I)+X*AP
!    5  AQ=Q3(I)+X*AQ
!       H=AP/AQ
!      ELSE
!       Y=1/X**2
!       H=(X-HF)*LOG(X)-X+C(4)+(C(1)+Y*(C(2)+Y*C(3)))/
!     1                                        ((C(5)+Y)*X)
!      ENDIF
!      DLGAMA=H
!      RETURN
!  101 FORMAT('NON-POSITIVE ARGUMENT  X = ',1P,E15.6)
!      END
        
      
FUNCTION DLGAMA(X)
  real(DP), intent(IN) :: X;
  CHARACTER*(*) NAME
  PARAMETER(NAME='ALGAMA/DLGAMA')
  real(DP) :: DLGAMA, H, Y, AP, AQ;
  real(DP) :: P1(7),Q1(7),P2(7),Q2(7),P3(7),Q3(7),C(5)
  real(DP), PARAMETER :: Z1 = D1, HF = D05, HF1 = 1._DP + HF;
  integer :: I;  
  CHARACTER*80 ERRTXT;  
  DATA P1 &
  /+3.842873656745991E+0_DP, +5.270689375300983E+1_DP, &
   +5.558404572351531E+1_DP, -2.151351357372570E+2_DP, &
   -2.458726172229242E+2_DP, -5.750089360304123E+1_DP, &
   -2.335909894951284E+0_DP/
  DATA Q1 &
  /+1.000000000000000E+0_DP, +3.373304790707074E+1_DP, &
   +1.938778403437713E+2_DP, +3.088295497342428E+2_DP, &
   +1.500683906489095E+2_DP, +2.010685134433395E+1_DP, &
   +4.571742028250299E-1_DP/
  DATA P2 &
  /+4.8740201396838636E+0_DP, +2.4884525168574076E+2_DP, &
   +2.1797366058895915E+3_DP, +3.7975124011525118E+3_DP, &
   -1.9778070769841646E+3_DP, -3.6929834005591282E+3_DP, &
   -5.6017773537803877E+2_DP/
  DATA Q2 &
  /+1.0000000000000000E+0_DP, +9.5099917418208938E+1_DP, &
   +1.5612045277928635E+3_DP, +7.2340087928948071E+3_DP, &
   +1.0459576594058959E+4_DP, +4.1699415153200231E+3_DP, &
   +2.7678583623804101E+2_DP/
  DATA P3 &
  /-6.880624009459425E+3_DP, -4.306996981957098E+5_DP, &
   -4.750459465343956E+6_DP, -2.942344593032234E+6_DP, &
   +3.632180493154257E+7_DP, -3.356778281454576E+6_DP, &
   -2.480436948828593E+7_DP/
  DATA Q3 &
  /+1.000000000000000E+0_DP, -1.421682983965146E+3_DP, &
   -1.555289028085353E+5_DP, -3.415251710801107E+6_DP, &
   -2.096962325580444E+7_DP, -3.454417509334395E+7_DP, &
   -9.160558286371317E+6_DP/
  DATA C &
  / 1.1224921356561E-1_DP,  7.9591692961204E-2_DP, &
   -1.7087794611020E-3_DP,  9.1893853320467E-1_DP, &
    1.3469905627879E+0_DP/

  !ENTRY DLOGAM(X)
  IF(X .LE. D0) THEN
    H = D0;
    WRITE(ERRTXT,101) X
   !CALL MTLPRT(NAME,'C304.1',ERRTXT)
    call dbg_assert_true(ERRTXT, X > D0);
  ELSE IF(X .EQ. D1 .OR. X .EQ. D2) THEN
    H = D0;
  ELSE IF(X .LE. HF) THEN
    Y = X + D1;
    AP = P1(1);
    AQ = Q1(1);
    DO I = 2,7
      AP = P1(I) + Y * AP;
      AQ = Q1(I) + Y * AQ;
    end do
    H = - LOG(X) + X * AP / AQ;
  ELSE IF(X .LE. HF1) THEN
    AP = P1(1);
    AQ = Q1(1);
    DO I = 2,7
      AP = P1(I) + X * AP;
      AQ = Q1(I) + X * AQ;
    end do
    H = (X-D1) * AP / AQ;
  ELSE IF(X .LE. D4) THEN
    AP = P2(1);
    AQ = Q2(1);
    DO I = 2,7
      AP = P2(I) + X * AP;
      AQ = Q2(I) + X * AQ;
    end do
    H = (X-D2) * AP / AQ
  ELSE IF(X .LE. D12) THEN
    AP = P3(1)
    AQ = Q3(1)
    DO I = 2,7
      AP = P3(I) + X * AP
      AQ = Q3(I) + X * AQ
    end do
    H = AP / AQ;
  ELSE
    Y = D1/(X * X);
    H = (X-HF)*LOG(X)-X+C(4)+(C(1)+Y*(C(2)+Y*C(3)))/ &
        ((C(5)+Y)*X);
  ENDIF
  DLGAMA = H;
  RETURN
101 FORMAT('NON-POSITIVE ARGUMENT  X = ',1P,E15.6)
  END function
end module